clear all; close all; clc; % clear workspace
load data.mat; % load data -- here variable name is 'npp'
% n.b. this code relies on tfn.m

%% fit to LSN

lnpp = log(npp); mu = mean(lnpp); sigma = std(lnpp); lnpp = (lnpp-mu)./sigma; % take logarithm, standardize, & record log-moments
gamma = skewness(lnpp); % standardized sample log-skewness, used for fitting initialization
gamma = max([-.99 gamma]); gamma = min([.99 gamma]); % so that initial guess is in range [-.99 .99]
delta = sign(gamma).*sqrt(pi/2*(abs(gamma).^(2/3))./(abs(gamma).^(2/3)+((4-pi)/2).^(2/3))); % define initial guess parameters
alpha = delta./sqrt(1-delta^2);
omega = std(lnpp)./sqrt(1-2.*delta.^2./pi);
xi = mean(lnpp)-(omega.*delta.*sqrt(2/pi));
[y_ecdf x_ecdf] = ecdf(lnpp); % define empirical CDF
y_ecdf = y_ecdf(2:end); x_ecdf = x_ecdf(2:end); % remove edgepoint ecdf function adds
lsncdf = normcdf((x_ecdf-xi)/omega)-2.*tfn((x_ecdf-xi)/omega,alpha); % define initial guess
l = max([lsncdf-y_ecdf; 0]); % calculate differences
u = max([y_ecdf-lsncdf; 0]);
ul = u+l;

%% iterate

n_it = 1000; % choose iteration number
snr = 25; % progressively smaller random jumps in parameter space
% can play around with different snr values

UL = [];
for i = 1:n_it; % randomly vary parameters and update if fit improves
    i./(3.*n_it)
    xi_it = awgn(xi,snr);
    omega_it = awgn(omega,snr);
    alpha_it = awgn(alpha,snr);
    lsncdf_it = normcdf((x_ecdf-xi_it)/omega_it)-2.*tfn((x_ecdf-xi_it)/omega_it,alpha_it);
    l_it = max([lsncdf_it-y_ecdf; 0]);
    u_it = max([y_ecdf-lsncdf_it; 0]);
    if l_it+u_it <= ul;
        xi = xi_it; omega = omega_it; alpha = alpha_it; ul = l_it+u_it; u = u_it; l = l_it;
    end
    UL(i) = ul;
end

snr = 50;
for i = 1:n_it;
    i./(3.*n_it)+1/3
    xi_it = awgn(xi,snr);
    omega_it = awgn(omega,snr);
    alpha_it = awgn(alpha,snr);
    lsncdf_it = normcdf((x_ecdf-xi_it)/omega_it)-2.*tfn((x_ecdf-xi_it)/omega_it,alpha_it);
    l_it = max([lsncdf_it-y_ecdf; 0]);
    u_it = max([y_ecdf-lsncdf_it; 0]);
    if l_it+u_it <= ul;
        xi = xi_it; omega = omega_it; alpha = alpha_it; ul = l_it+u_it; u = u_it; l = l_it;
    end
    UL(n_it+i) = ul;
end

snr = 100;
for i = 1:n_it;
    i./(3.*n_it)+2/3
    xi_it = awgn(xi,snr);
    omega_it = awgn(omega,snr);
    alpha_it = awgn(alpha,snr);
    lsncdf_it = normcdf((x_ecdf-xi_it)/omega_it)-2.*tfn((x_ecdf-xi_it)/omega_it,alpha_it);
    l_it = max([lsncdf_it-y_ecdf; 0]);
    u_it = max([y_ecdf-lsncdf_it; 0]);
    if l_it+u_it <= ul;
        xi = xi_it; omega = omega_it; alpha = alpha_it; ul = l_it+u_it; u = u_it; l = l_it;
    end
    UL(2.*n_it+i) = ul;
end
lsncdf = normcdf((x_ecdf-xi)/omega)-2.*tfn((x_ecdf-xi)/omega,alpha);

%%

XI = xi; % final xi value
OMEGA = omega; % final omega value
ALPHA = alpha; % final alpha value
V = u+l; % final V value

clearvars -EXCEPT XI OMEGA ALPHA V;